Traveling Dark Solitons in Superfluid Fermi Gases 
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Families of dark solitons exist in superfluid Fermi gases. The energy-velocity dispersion and number of 
depleted particles completely determines the dynamics of dark solitons on a slowly-varying background density. 
For the unitary Fermi gas we determine these relations from general scaling arguments and conservation of local 
particle number. We find solitons to oscillate sinusoidally at the trap frequency reduced by a factor of 1/ \/3. 
Numerical integration of the time-dependent Bogoliubov-de Gennes equation determines spatial profiles and 
soliton dispersion relations across the BEC-BCS crossover and proves consistent with the scaling relations at 
unitarity. 
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Dark solitons are elementary nonlinear excitations that play 
a key role in understanding complex dynamics of superflu- 
ids |Q], 12]. Superfluid Fermi gases have only recently become 
accessible experimentally and their nonlinear wave dynamics 
are largely unexplored yL|4(]. These systems offer the intrigu- 
ing possibility to tune between the perturbatively accessible 
regimes of Bose-Einstein condensation (BEC) of preformed 
pairs and Bardeen-Cooper-Schrieffer (BCS) superfluidity and 
a strongly correlated regime of unitarity-limited interactions. 
While the existence and properties of dark solitons in the BEC 
regime can be inferred from the solutions of Gross-Pitaevskii 
(GP) mean-field theory and experiments with atomic BECs, it 
is an outstanding question what happens outside this regime. 
So far, only numerical solutions for stationary dark solitons 
within Bogoliubov-de Gennes (BdG) mean-field theory have 
been available [5]. 

In this work, we report theoretical results supporting the 
existence and detailing the properties of a family of travel- 
ing (grey) solitons that are parameterized by their velocity of 
propagation v s . We are aware of parallel efforts to understand 
soliton dynamics in trapped Fermi gases [6] and to determine 
grey soliton profiles [7]. For the unitary gas, we find a closed 
analytic form of the energy-velocity dispersion relation that is 
fully determined from a set of general assumptions: (a) Upon 
adiabatic change of the environment, the soliton can adjust its 
dynamical state and consistently conserve locally both energy 
and particle number, (b) Energy and particle number vanish 
as v s approaches the speed of sound, (c) The superfluid order 
parameter has a well defined phase step across the soliton that 
also vanishes under the conditions of (b). 

Assumption (a) means that solitons can move adiabatically 
between regions of different background density without dis- 
integration and radiation. This non-trivial property is known 
to be true for GP solitons (§]■ For the unitary gas the as- 
sumptions are supported by the excellent agreement found be- 
tween the analytic dispersion and our numerical results based 
on BdG mean-field theory (see Figs.QJ) and|2}3,d). The mean- 
field calculations further allow us to obtain spatial soliton pro- 
files and dispersion relations for arbitrary interactions outside 
the unitarity limit. 
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FIG. 1. (Color online) (a) Soliton energy E s as a function of velocity 
v s for different couplings: (a) r\ — —0.5 (BCS), (b) 77 = (uni- 
tary), and (c) r\ = 1 (BEC), where the solid lines show the analytical 
relations 10 and E < f p — ^hnsc [l — u 2 /c 2 ] 3 ^ 2 , respectively. 



As a dynamical consequence we are able to predict os- 
cillations of dark solitons in a harmonically trapped Fermi 
gas. While for BECs, the oscillation frequency was predicted 
OSl] and observed Jiol [Till to be reduced from the trapping 
frequency uj t by a factor of l/\/2 « 0.707, we find the 
oscillation frequency further reduced across the BEC-BCS 
crossover. The BdG calculation yields w/w t =0.480, 0.572, 
and 0.687 for 77 = -0.5 (BCS regime), (unitary) and 1 (BEC 
regime), respectively, where 77 = l/(kFd), the Fermi wave 
number kp = {i^n) 1 ^ parametrizes the density n, and a is 
the s-wave scattering length. Our analytic theory for the uni- 
tary case of 77 = predicts uj/ut = l/v3 ~ 0.577, which is 
in excellent agreement with the numerical data. 

Let us consider a superfluid Fermi gas with a soliton that 
is localized along the z direction in a region small compared 
to the system length L on a homogeneous background. We 
can extract the scaling of the system energy with density us- 
ing the inverse Fermi wave number hp 1 and the Fermi en- 
ergy Ep = h 2 k F /(2m) as units of length and energy, re- 
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FIG. 2. (Color online) Soliton properties at unitarity r\ — as a 
function of velocity: (a) the magnitude of the order parameter in the 
center, (b) the phase difference 8(f) = arg[A(oo)] — arg[A(— oo)], 
(c) the density in the center, and (d) the number of particles in the 
soliton. Dots show numerical data and lines plot Eqs. ® ar, d ©, 
respectively. 



spectively. We may express the grand canonical energy of 
the superfluid Fermi system E 1 = (H — fiN) = [ehkpL + 
£ + (D((kpL)~ 1 )]k F AEp, where A is the transverse area, Eh 
the dimensionless background energy and E s = £k 2 F AEp 
is identified as the soliton energy. Physically, the soliton en- 
ergy E s (p,,v s ,a) depends on three independent parameters, 
where v s is the propagation velocity. Dimensionless £ , how- 
ever, may only depend on the two dimensionless parameters of 
velocity v = v s /vp and coupling strength r/. The dependence 
on the chemical potential p, is implicit through the density, 
which determines kp. 

A special case arises for unitarity, where a —> oo and 77 = 
becomes independent of /1. Anticipating that E s is an even 
function of the velocity, we consider £(v 2 ) as a function of 
v 2 . Employing the equation of state of the unitary Fermi gas 
fi = (1 + j3)Ep, where (3 is the many-body parameter, we can 
extract the dependence on fi 



E s (fx,v s ) = fi 2 B£(v 2 ) 



(1) 



where d 2 =v%(l + /3)m/(2/x) and B = 2Am/[(l + f3)H} 2 . 
As a localized wave form, a soliton is characterized not only 
by its energy but also by its particle number N s = J (n s — 
no)d 3 r — —dE s /dfi, where n s is the soliton density and no 
is the background density [8]. We find 



N s = v 2 -{l+P)B£\i 2 ) - 2nB£(v% 



(2) 



where £' — d£/dv 2 . The particle number and energy are thus 
completely determined by the same function £ (v 2 ). We are 
now going to determine this function from the assumptions 
(a) - (c). 

Assumption (a) requires both E s and N s to be constants of 
the motion through inhomogeneous density. Since we are con- 



sidering purely one-dimensional motion of the soliton, treat- 
ing it as a quasiparticle, there is at most only a single inde- 
pendent constant of the motion. The condition that N s and E s 
have identical contours in phase space leads to the condition 



dN, dE, ON, dE q 



d/j, dv s dv s djj. 



(3) 



where \x and v s represent the quasiparticle's coordinate and 
momentum, respectively. It follows from Eqs. (Q]) and (0 that 
the third derivative of £ vanishes identically and that £ (v 2 ) 
can be parameterized by 



£(i 2 ) =e(v 2 



(4) 



where e and v are yet undetermined parameters. The func- 
tional form © already has important implications for soliton 
oscillations in a trapped gas: 

Requiring dE s /dt = 0, we find the Newtonian equation of 
motion for the soliton 



d/j, 
dz 



m(l + j8) 



0. 



(5) 



In the case of harmonic trapping and under validity of the 
Thomas Fermi approximation, we can write fi(z) = [iq — 
muj 2 z 2 /2 and Eq. (O reduces to a harmonic oscillator. The 
frequency u)/u)t = v/\/l + (3 is independent of amplitude! 

The parameter v can be determined from assumption (b): 
From Eq. (0]l we find that the energy and particle number 
vanish when the dimensionless velocity v reaches the criti- 
cal value v. We expect this to happen at the speed of sound, 
which takes the value c = v/(l + j3) /3vp. This leads to 
V = y(l + /3)/3 and yields the oscillation frequency uj /uj t = 
1/a/3. We have thus derived the oscillation frequency of a 
dark soliton in a harmonically trapped unitary gas from the 
assumptions (a) and (b). 

The remaining coefficient e can be determined from the re- 
lation between the physical momentum of the soliton p s — 
mN s v s and the canonical momentum p c , which is defined by 
dEg/dpclft — v s . The difference between the two quantities 
accounts for the counterfiow that would have to occur in a 
toroidal system to compensate for the phase difference 5(f) in 
the superfluid order parameter across the soliton [12]. For the 
superfluid Fermi gas the counterfiow term was recently found 
by Pitaevskii |0] 



Ps~Pc = tin^n - 5<f>)/2, 



(6) 



where ni = nA = k F A/ (3ir 2 ) is the one-dimensional den- 
sity. From Eqs. (Q~|) and (01 we evaluate the difference using 

p c — J v~ l dE s /dv s dv s = —hni6ir 2 ev[v 2 — w 2 /3] to yield 



Ps — Pc = fvni^K 2 ev 2 v. 



(7) 



Comparing Eqs. (|7]i and ©, we find that the phase differ- 
ence varies linearly with velocity in contrast to the GP soli- 
ton, where cos(<50 GP /2) = v s /c GP . Fixing the remaining 
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constant e by requiring the phase step to vanish at the speed 
of sound [assumption (c)], we find 



— , 6<p = 7r(l-v s /c). 
an 



(8) 



Thus, the energy and particle number dispersion (shown as 
full lines in Figs, lb and 2d, repsectively) as well as the phase 
step for the family of dark solitons in the unitary gas (shown 
in Fig. 2b) are obtained without any free parameters. The 
success of this derivation shows, that the assumption (a) of 
particle-number conservation under quasiparticle motion is 
consistent with the universal scaling relations of the unitary 
Fermi gas. 

We have not yet proven that dark solitons exist. Within 
the realm of mean-field theory, this can be done by finding 
self-consistent solutions of the BdG equations. In addition to 
testing the stated assumptions against a physical theory, this 
allows us to determine spatial profiles as well as dispersion 
relations outside the unitary regime. 

We now more generally consider a Fermi gas with equal 
density for two spin components at zero temperature. The 
time-dependent BdG equations provide a convenient mean- 
field theory of the BEC-BCS crossover fl 

h A(r,ty 
,A*(r,i) -h 



ihdf 



u v {r,t) 
v y {r,t) 



u v {r,i) 
v u (r,t) 



(9) 



where h = |-V 2 - 

2m 



/j and u and v are space- 
and time-dependent quasi-particle amplitudes satsifying 
/ d 3 r [u*(r,t)u„/(r,t) + v*(r,t)v u *(r,t)} = § vv >. The prob- 
lem simplifies to a time-independent eigenvalue problem 
when we seek soliton solutions of the superfluid order pa- 
rameter of the form A(z,t) — A(z — v s t) — A(£) and 



write v v (r, t) — {LA) 



iE, 



■"' 7 %,n(0 and 



likewise for u. The energies E P:1l are the eigenvalues of 
the resulting time-independent BdG equation, which contain 
the soliton velocity v s as a parameter. The transverse mo- 
mentum p is discretized according to the transverse area A 
of the computational box. The above equations must be 
solved together with the equation for the order parameter 
A (£) = -sEp.nVKl^K) in a self-consistent way. 
The density is then given by = 2 ^ p n |i>p, ra (£)| 2 . All 
sums are restricted to < E p>n < E c where E c is a high 
energy cutoff. The coupling strength g relates to the scat- 
tering length a through the cutoff-dependent renormalization 
1/g = m/{Airh 2 a) - l/2e„ 0H where e„ is the 

energy of the state v in the normal phase. Open boundary con- 
ditions are imposed by ensuring the hermiticity of the matrix 
and the proper symmetry of A(£) in a self-consistent way. We 
have implemented a generalized secant (Broyden's) method to 
find self-consistent solutions (with very small self-consistency 
error of 10 -8 or better) from an initial profile with nontrivial 
phase structure. 

The spatial structure of the dark soliton solutions is shown 
in Figs. [3] and |4] for three values of the interaction parameter 
rj. The main feature is a density notch and dip in the order pa- 
rameter, which become shallower with increasing velocity for 
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FIG. 3. (Color online) The spatial structure of the soliton order pa- 
rameter A at different velocities by its imaginary part (a) and magni- 
tude (b) in the BCS regime at 77 = —0.5 (la,b), unitarity limit 77 = 
(2a,b), and BEC regime rj = 1 (3a,b). 



all interacting regimes. In the BEC regime, the profiles most 
closely resemble the GP dark solitons with constant imagi- 
nary part (with appropriate choice of a global phase) and bell- 
shaped, const — sech 2 (z/l), density notch. This is expected as 
the BdG equation converges toward the GP equation for large 
positive 77 ni5ll . The imaginary part of the order parameter 
clearly develops structure for finite velocities in the unitarity 
and BCS regime, which is a striking new feature of the BdG 
grey solitons as seen in panels (la) and (2a) in Fig. [3] 

While the length scale for the BEC soliton is £ = 
h/ {m^J c 2 — vf) from GP theory, there is no clear evidence 
in our data for a velocity dependence of the length scale in 
the unitarity limit. There we expect on general grounds that 
the only length scale is kp 1 . In the BCS regime, we expect 
small scale Friedel oscillations with size 2/c^ 1 , and a second 



length scale in the Cooper pair size £, 



evaluates to about hkp 1 for 77 = —0.5 |5 



= Hvf/Aq, which 
. From Fig.fSJ panel 



(la) and our experience with boundary effects, it appears that 
there is an additional velocity dependence and the total size 
grows with increasing velocity. 

Relevant velocity scales for the problem are the speed of 
sound c = n(d n / dn) / m and the pair breaking velocity 
/i. We consistently found it difficult 



A 2 
^0 



to converge to self-consistent solutions approaching these ve- 
locities from below and thus assume that soliton solutions ex- 
ist only below v c = min(c, v sp ), which is also the critical 
velocity for dissipationless motion of infinitesimal impurities 
J3L rT^ll . v c takes a maximum around unitarity, where c s» v sp . 
For 77 > (BEC to unitarity), c < v sp and thus solitons are 
limited by the speed of sound. In the BCS regime, pair break- 
ing dominates and v sp limits soliton propagation. 

The dimensionless soliton energy as a function of the ve- 
locity is shown in Fig. [T] In all three regimes the energy is 
positive with negative curvature, which supports the under- 
standing of a dark soliton as a quasiparticle with negative ef- 
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FIG. 4. (Color online) The density profile at different velocities for 
(a) rj = -0.5 (BCS), (b) r/ = (unitarity) and (c) jj = 1 (BEC). 



fective mass. In the unitarity limit, the numerical data fits 
beautifully with the analytical result (0|, where we have used 
1 + [$ = 0.6082 consistent with our BdG calculation of the 
chemical potential. The GP formula for the energy shown in 
panel (c) is strictly valid only in the limit of large rj, which 
explains the small deviations from the numerical data. 

Additional properties of dark solitons are shown in Fig. [2] 
for the unitarity limit. The order parameter at the soliton cen- 
ter shown in panel (a) intriguingly appears to vary linearly 
with velocity, as would be the case for GP solitons. The phase 
step and particle number agree very well with the analytical 
predictions based on Eqs. (0| and ©. We have further tested 
the assumption (a) that N s is a constant of the motion by di- 
rectly checking Eq. (01, which is fulfilled within our expected 
numerical errors. 

Finally, we consider small amplitude soliton oscillations in 
a harmonically trapped Fermi gas beyond the unitarity regime. 
From the Thomas-Fermi and local density approximation we 
have E s (v s , //rj — ™ t 2 z s 2 /2) = const and taking a time deriva- 
tive obtain 



1 8E t 

v s dv s 



z s + N s muj t z. 



0. 



(10) 



Noting that both E s and N s should be even functions of ve- 
locity, Eq. ( TTOb describes a harmonic oscillator. Defining the 
effective mass M s = vJ 1 ^ s -\ Vs -o, we find the frequency of 
small oscillations 



'miV s (0) 



(11) 



mass mN s and effective mass M s . Both N s (0) and M s are 
negative and can be easily extracted from our numerical data. 
We find the oscillation frequencies u)/ui t =0.480, 0.572, and 
0.687 for r/=-0.5, and 1 respectively. 

The oscillation frequency of solitons thus decreases signif- 
icantly from the BEC towards the BCS regime in agreement 
with time-dependent simulations J6J]. In the unitarity regime 
we were able to determine the dispersion relation in closed 
form starting from a small number of global assumptions. In 
particular we find that the oscillation frequency does not de- 
pend on the many-body parameter f3 in the equation of state 
or other details of the soliton solutions. Dark solitons thus of- 
fer the opportunity to explore complimentary properties of the 
unitary gas to previous studies that pinpointed the equation of 
state [31]. The important questions of stability of dark soli- 
tons against strong quantum fluctuations and the consistency 
of local energy and particle number conservation deserve to be 
studied beyond mean-field theory and thus make an excellent 
subject for future experimental investigation. 
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